function [Ea,Eu,AddEa,AddEu] = exogenize(S,M,Yxx,Ea,Eu)
% exogenize  [Not a public function] Compute add-factors to endogenised shocks.
%
% Backed IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

%--------------------------------------------------------------------------

if isempty(S.Anch) || ~any(S.Anch(:))
    return
end

ny = size(S.Z,1);
nxx = size(S.T,1);
ne = size(Ea,1);
nPer = size(Ea,2);

yxxTune = S.Tune(1:ny+nxx,:);
yxxAnch = S.Anch(1:ny+nxx,:);
eaAnch = S.Anch(ny+nxx+(1:ne),:);
euAnch = S.Anch(ny+nxx+ne+(1:ne),:);
eaWght = S.Wght(1:ne,:);
euWght = S.Wght(ne+(1:ne),:);

% Last exogenized and endogenized data points.
% findLastFn = @(X) max([0,find(any(any(X,3),1),1,'last')]);
% lastExg = findLastFn( yxxAnch(1:ny+nxx,:) );
[lastEndgA,lastEndgU] = utils.findlast(eaAnch,euAnch);

% Compute prediction errors.
% pe : = [ype(1);xpe(1);ype(2);xpe(2);...].
pe = yxxTune(yxxAnch) - Yxx(yxxAnch);

% Compute add-factors that need to be added to the current shocks.
if size(M,1)==size(M,2)
    
    % Exactly determined system
    %---------------------------
    upd = M \ pe;

else
    
    % Underdetermined system (larger number of shocks)
    %--------------------------------------------------
    d = [ ...
        eaWght(eaAnch); ...
        euWght(euAnch) ...
        ].^2;
    nd = length(d);
    P = spdiags(d,0,nd,nd);
    upd = simulate.linear.updatemean(M,P,pe);
    
end

nnzEa = nnz(eaAnch(:,1:lastEndgA));
ixEa = eaAnch(:,1:lastEndgA);
ixEu = euAnch(:,1:lastEndgU);

if issparse(Ea)
    [row,col] = find(ixEa);
    AddEa = sparse(row,col,upd(1:nnzEa),ne,nPer);
    Ea = Ea + AddEa;
else
    AddEa = zeros(ne,lastEndgA);
    AddEa(ixEa) = upd(1:nnzEa);
    Ea(:,1:lastEndgA) = Ea(:,1:lastEndgA) + AddEa;
end

AddEu = zeros(ne,lastEndgU);
AddEu(ixEu) = upd(nnzEa+1:end);
Eu(:,1:lastEndgU) = Eu(:,1:lastEndgU) + AddEu;

end
